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Abstract 

Background: Ticks are the most important pathogen vectors in Europe. They are known to be influenced by 
environmental factors, but these links are usually studied at specific temporal or spatial scales. Focusing on Ixodes 
ricinus in Belgium, we attempt to bridge the gap between current "single-sided" studies that focus on temporal or 
spatial variation only. Here, spatial and temporal patterns of ticks are modelled together. 

Methods: A multi-level analysis of the Ixodes ricinus patterns in Belgium was performed. Joint effects of weather, 
habitat quality and hunting on field sampled tick abundance were examined at two levels, namely, sampling level, 
which is associated with temporal dynamics, and site level, which is related to spatial dynamics. Independent 
variables were collected from standard weather station records, game management data and remote sensing-based 
land cover data. 

Results: At sampling level, only a marginally significant effect of daily relative humidity and temperature on the 
abundance of questing nymphs was identified. Average wind speed of seven days prior to the sampling day was 
found important to both questing nymphs and adults. At site level, a group of landscape-level forest fragmentation 
indices were highlighted for both questing nymph and adult abundance, including the nearest-neighbour distance, 
the shape and the aggregation level of forest patches. No cross-level effects or spatial autocorrelation were found. 

Conclusions: Nymphal and adult ticks responded differently to environmental variables at different spatial and 
temporal scales. Our results can advise spatio-temporal extents of environment data collection for continuing 
empirical investigations and potential parameters for biological tick models. 
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Background 

Ixodes ricinus is a hard tick species (Acari: Ixodidae) 
widely distributed in Europe. It is capable of transmitting 
a number of pathogens to humans and livestock. For ex- 
ample, Lyme borreliosis, caused by the bacterium Borre- 
lia burgdorferi, is the most frequent tick-borne disease 
of humans in temperate zones [1]. Louping ill, an acute 
viral disease primarily of sheep, can result in a mortality 
of up to 60% for infected sheep. Understanding the 
spatio-temporal dynamics of I. ricinus can be of public 
health and economic importance. 
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The current literature on tick patterns documents two 
types of associations: (i) between the extrinsic climatic 
factors and temporal dynamics and (ii) between the 
physical/biotic environmental factors and spatial dynam- 
ics. Studies on the temporal dynamics of ticks generally 
investigate tick phenology, which refers to the seasonal 
cycle of development and activity of the life stages of 
tick and the impacts from seasonal and inter-annual cli- 
mate variations. Tick survival and development rates 
can be influenced by relative humidity and temperature 
[2,3] and the host-seeking activity can be affected by day 
length and vapour pressure deficit [4]. Kurtenbach et al. 
[5] proposed that the phenology of ticks shapes the tem- 
poral patterns of tick-borne diseases. Studies on the 
spatial dynamics of ticks often focus on the suitability of 
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habitat. Preferential habitats of ticks and many of their 
hosts occupy a large fraction of rural landscapes. Habitat 
types, regarding the landscape composition, landscape 
configuration, soil conditions, bedrock geology, topog- 
raphy, etc., have been largely and empirically tested for 
their associations to the dynamics of ticks. Forests are 
the most suitable habitats for ticks in most studies in 
North America (I. scapularis) and Europe (/. ricinus) [6]. 
Their configurations have been identified as key deter- 
minants of the spatial dynamics of ticks [7,8]. The abun- 
dance of questing (I. scapularis) ticks is positively 
associated with sandy or sandy loam soils and sediment- 
ary bedrock [9]. In the Czech Republic [10], Switzerland 
[11] and Spain [12], topographic variables also poten- 
tially influence the local tick abundance. Landscape can 
shape the patterns of host dynamics. Key hosts such as 
bank voles and deer inhabit agrarian and forested land- 
scape [5]. Moreover, ecotones between forest and culti- 
vated fields, pastures and other grassland areas can be of 
high value to ticks and some host species [13]. The rela- 
tionship between tick abundance and landscape struc- 
ture, however, has been sparsely investigated. 

However, current investigations suffer from several 
limitations. Firstly, researchers have rarely investigated 
tick patterns both temporally and spatially. One study 
examined the influence of weather and biotic (host) fac- 
tors in northern Italy [14], but did not take physical en- 
vironmental factors into account. Secondly, there is 
limited understanding of the influence of forest fragmen- 
tation on the spatial dynamics of ticks. A recent regional 
level study demonstrated a positive effect of forest edge 
on questing tick abundance in northern Belgium [15]. 
Forest patch size was also found negatively associated to 
(/. scapularis) tick abundance in southeastern New York 
[8]. In northern Spain, high tick abundance was found in 
stepping-stone' patches that cause significant changes in 
connectivity when removed from the landscape [7]. 
Other than these patch-level effects, a number of 
landscape-level characteristics describing the spatial con- 
figuration of patches (e.g. neighbouring distance, aggre- 
gation level, fractal dimension) can be also important as 
they may shape the local host exchange patterns. Thirdly, 
the influence of wildlife management needs to be further 
examined. Effects of deer exclusion by fencing have been 
empirically studied in Scotland [16] and Italy [17]. Hunt- 
ing, however, has received less attention thus far. Hunt- 
ing can directly influence the local populations of large 
vertebrates (e.g. deer), which are key reproduction hosts 
for adult female ticks. It can be important to explore 
whether hunting can be used as a strategy to reduce the 
local tick abundance, provided the spatial dataset of wild- 
life are sometimes too difficult to get. 

To address these limitations, a multi-level analysis was 
carried out to study the spatio-temporal dynamics of /. 



ricinus abundance in Belgium. A multi-level approach is 
considered suitable for the present study because of the 
clustered nature of our data, i.e. repeated tick samplings 
in most individual sites. Despite a growing interest in 
this approach in the field of vector-borne diseases [18], 
it has, to our knowledge, never been used in previous 
tick-related studies. In the present study, joint effects of 
weather, habitat quality and wildlife management factors 
were examined at two levels: (i) sampling level, which is 
associated with temporal tick dynamics, and (ii) site 
level, which is related to spatial tick dynamics. The main 
objective was to determine in each level the key factors 
and to examine possible cross-level effects to the spatio- 
temporal dynamics of questing ticks in Belgium. 

Methods 

Study sites 

From 2009 to 2010, 125 samplings were carried out in 
51 sites throughout Belgium (Figure 1): 36 sites were vis- 
ited in 2009 and 27 sites were visited in 2010, among 
which 12 sites were visited in both years. Sites were ini- 
tially selected based on a 25 km * 25 km regular grid 
across Belgium, from which unsuitable habitats (such as 
built-up areas) were excluded. In the field, some sites 
were moved, in order to be placed in a suitable forest 
stand. All sampling sites were in broad-leaved or mixed 
forest stands. Most of these forest stands are adjacent to 
cropland and pasture only, while a few of them were also 
adjacent to urban and water areas. 

Ixodes ricinus sampling scheme 

Standard blanket-dragging was adopted to estimate the 
abundance of the questing ticks from March to Novem- 
ber in 2009 and from March to September in 2010. 
These timings were selected because the 'tick season' 
generally lasts from March to October in Belgium [19]. 
During each sampling, aim 2 piece of white flannel was 
pulled over the short grass or litter along a 500 m tran- 
sect. High ground vegetation such as purple moor grass 
and bracken fern was avoided. The blanket was exam- 
ined every 20 m and nymphs and adults were preserved 
in alcohol for future laboratory identification. The iden- 
tification was done on the basis of morphologic features 
with a stereo-microscope by using standard taxonomic 
keys [20]. All samplings were performed between 9:30 
and 13:00 and rainy days were excluded. At the begin- 
ning of each sampling, coordinates were recorded by 
military GPS (spatial resolution 1 m). 

Environmental variables 

Environmental variables were collected at two levels: (i) 
at the level of sampling event ("sampling level"), focusing 
on the influence of seasonal weather variations on the 
temporal patterns of ticks; (ii) at site level, considering 
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the role of landscape (composition and configuration) 
and soil conditions on the quality of tick habitat, as well 
as hunting activity that affects the population of host 
animals. 

Weather data was acquired from the Royal Institute of 
Meteorology for all working stations covering the period 
of sampling. Measures were available from 50 meteoro- 
logical stations for daily maximum and minimum tem- 
peratures and precipitation and from 25 stations for 
daily wind speed and relative humidity. The average dis- 
tance between sampling sites and the closest station for 
temperature and precipitation is 17.3 km and that for 
wind speed and relative humidity is 12.3 km. For each 
sampling session, we selected the nearest meteorological 
stations for daily maximum and minimum temperatures, 
precipitation, wind speed and relative humidity. In 
addition, average daily vapour pressure deficits (VPD) 
were estimated from observations of temperature and 
relative humidity, following [21]: VDP = es(Ta) x (1-RH), 
where es(Ta) is the saturation vapour pressure function of 
mean temperature Ta and RH is the relative humidity. 
Mean values of these variables in 3 and 7 days prior to the 
sampling day were also calculated. 

A land cover classification was generated based on 
Landsat satellite images (combined from two captures 
on 2000/09/11 and 2001/05/23; resolution 29.18 m). 
Landscape composition was measured in 500 m-, 
1000 m- and 1500 m- buffers around each sampling site. 
The radius of 500 m was chosen because it was the total 
length for one field sampling. Radii of 1000 m and 



1500 m were selected to examine the landscape effects 
at larger extents. The proportions of potential habitat 
types (i.e. broad-leaved forest, coniferous forest, heath- 
land and agricultural areas of crops and pastures) were 
calculated. The configuration of forests was examined 
within the same buffers using several indices [22]: the 
number of forest patches, edge density (the sum of the 
lengths of all edge segments involving the forest type 
divided by the total landscape area), aggregation index 
(aggregation levels of forest patterns), Euclidean mean 
nearest-neighbour distance (mean distance to the near- 
est neighbouring forest patches), area- weighted mean 
patch fractal dimension (overall complexity and frag- 
mentation of forest patches) and area-weighted mean 
shape index (overall shape complexity of forest patches). 

Two digital soil maps, for Wallonia (by "Unite Sol- 
Ecologie-Territoire, FUSAGx") and for Flanders (by 
"Vlaamse Landmaatschappij"), were acquired. Informa- 
tion on soil texture and drainage was extracted from 
both maps. For each sampling site, soil textures were 
coded as percentages of sand, silt, clay and gravel. Drain- 
age conditions were indexed from 0 to 8 from excessive 
to poor drainage across Belgium, while values of only 0 
to 6 were observed for the sampling sites. 

Data on roe deer shot in 2009 were collected. In 
Wallonia, data was provided for each forest-management 
administrative unit (i.e. "cantonnement") by "Direction 
Generale des Ressources Naturelles et de TEnvironne- 
ment". In Flanders, data was provided by "Instituut voor 
Natuur- en Bosonderzoek" and recorded in each 
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administrative unit (i.e. "wildbeheereenheden", a unit 
that is much smaller than "cantonnement"). For each 
administrative unit ("cantonnement" for Wallonia or 
"wildbeheereenheden" for Flanders), numbers of roe 
deer shot were divided by the forested area to get the 
number of shot roe deer per km 2 of forest. The data 
may be biased as officials from the two regions focused 
on different spatial scales and may put different weights 
on information sources when making estimations. To 
weaken the effects of such uncertainty, we classified the 
number of roe deer shot per km 2 of forest into five 
categories using an interval of 50. 

Statistical analysis 

The abundances of nymphs and adults were analysed 
separately using multi-level approaches. All statistical 
analyses were conducted with R (version 2.14.1) [23]. 

In the first step, a test for overdispersion (i.e. likeli- 
hood ratio of a fitted Poisson model against a fitted 
negative binomial model) indicated that the data was 
overdispersed and that using a negative binomial distri- 
bution was necessary. A non-parametric Kruskal-Wallis 
one-way analysis of variance was used to test the neces- 
sity to use random intercept models. The between-site 
variance of intercepts was significantly different from 
zero (nymphs: chi-squared = 89.3, df = 50, P < 0.001; 
adults: chi-squared = 83.5, df = 50, P < 0.005). Thus, 
models that allow for random intercepts are better sui- 
ted than those with fixed intercepts. 

In the second step, single level regressions (i.e. bivari- 
ate negative binomial regressions) were performed to se- 
lect significant variables (P < 0.1) for each level. At 
sampling level, the dependent variable was the tick sam- 
pling result per site and the weather data were examined 
as independent variables. At site level, the dependent 
variable was defined as the mean value of sampled ticks 
in each site. Landscape composition and configuration, 
soil and hunting data were analysed as independent 
variables. 

In the third step, a set of multi-level negative binomial 
models with random intercepts were developed for 
model selection. Instead of evaluating a single model, a 
model selection approach offers a way to draw infer- 
ences from competing models by evaluating their rela- 
tive support in the data [24]. An Automatic 
Differentiation (AD) model builder approach in R, pack- 
age glmmADMB [25], was used to build the multi-level 
models. We generated nine groups of models by adding 
up to two variables for each level. For each group, mod- 
els were built in an automated fashion by looping all 
combinations of variables selected in the second step. 
Models were rejected if not all predictors were signifi- 
cant at <0.1 or if significant colinearity was observed 
(Pearson correlation coefficients between independent 



variables >0.7 and P < 0.1) among variables. We sorted 
the remaining models by their improvements on the 
model fit in comparison to intercept-only models as 
assessed by the change in Akaikes Information Criterion 
(AIC) (AAIC = AIC of intercept-only model-AIC of 
model with variables). Variables were ranked according 
to their relative importance (by summing the Akaike 
weights across the models that include the variable [26]) 
in the above-average models (i.e. models for which 
AAICs > averaged AAIC value of the remaining mod- 
els). The residuals of the top five models with the great- 
est model fitness improvement (the largest AAICs) were 
checked for spatial autocorrelation in model residuals, 
the significance of random slope and cross-level effects. 

Results 

Spatio-temporal variation of questing nymphs and adults 

A total of 22 400 I. ricinus ticks were collected, includ- 
ing 20 725 nymphs and 1 675 adults (828 males and 847 
females). More samplings were conducted in Flanders 
(93 times) than in Wallonia (32 times). Detailed seasonal 
patterns of sampling for each province are presented in 
Table 1. In general, I. ricinus ticks were found in all ten 
provinces in Belgium. Tick abundance per sampling was 
calculated for all ten provinces: in Flanders the sampled 
tick abundance was higher in East Flanders, Flemish 
Brabant and Limburg; while in Wallonia, that was higher 
in Namur, Hainaut and Luxembourg. 

The temporal variation of ticks in 2009 and 2010 is 
presented in Figure 2. In both years, the abundance of 
nymphs started to increase in spring, then fluctuated 
and reached a first peak in summer (July for 2009 and 
June for 2010)-one month before the hottest month. 
The second peak came in late summer and early au- 
tumn-one month after the hottest month. The associ- 
ation between the sampling results and meteorological 
conditions appears to be varying throughout the two 
sampling years. In both springs, questing tick abundance 
increased while temperature increased. The decline of 
the nymph abundance in the summer of 2009 can be 
associated to a peak in vapour pressure deficit or to a 
drop in maximal relative humidity. In the summer of 
2010, nymph abundance peaked while vapour pressure 
deficit peaked and relative humidity dropped. Finally, in 
both autumns, questing tick abundance decreased while 
temperature and vapour pressure deficit decreased but 
relative humidity increased. 

Single level regressions 

Results for regressions at sampling level and for regres- 
sions at site level are presented in Tables 2 and 3 together 
with descriptive statistics of the collected variables. At 
sampling level, abundance of nymphs was found signifi- 
cantly associated with maximum temperature, relative 
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Table 1 Ixodes ricinus sampled between 2009 and 2010 from ten provinces in Belgium (S = number of sampling; N.N = 



number of nymphs; N.A = 


number of adults; 


//_// 


not sampled) 










Province 


Year (20-) Spring (Mar-May) 
S N.N N.A 


Summer (Jun-Aug) 
S N.N N.A 


Autumn (Sep-Nov) 
S N.N N.A 


Nymphs per 
100 m 2 


Adults per 
100 m 2 


Flanders Antwerp 


09 


4 


627 


103 


1 


245 


15 


3 194 14 


26.7 


3.3 




10 


3 


513 


22 


4 


578 


82 


2 217 16 


29.1 


2.7 


East Flanders 


09 


4 


219 


84 


2 


427 


50 


_ 


21.5 


4.5 




10 


5 


1683 


47 


2 


541 


32 


_ _ 


63.5 


2.3 


Flemish Brabant 


09 


10 


2769 


278 


6 


2393 


109 


11 1458 204 


49.0 


4.4 




10 


8 


2420 


161 


6 


1802 


44 


2 529 30 


59.4 


2.9 


Limburg 


09 


2 


153 


16 


1 


298 


45 


2 334 37 


31.4 


3.9 




10 


1 


455 


14 


3 


818 


33 


1 142 7 


56.6 


2.2 


West Flanders 


09 


2 


0 


1 


2 


23 


8 




1.2 


0.5 




10 


3 


2 


1 


3 


3 


0 


_ _ 


0.2 


0.0 


Wallonia Hainaut 


09 








1 


125 


3 


_ _ 


25.0 


0.6 




10 


3 


65 


11 


2 


5 


1 


- 


2.8 


0.5 


Liege 


09 


1 


10 


0 


3 


144 


8 




7.7 


0.4 




10 


3 


140 


9 


1 


26 


1 




8.3 


0.5 


Luxembourg 


09 
10 


1 
1 


46 
86 


10 

3 


3 


196 


21 


_ 


12.1 
17.2 


1.6 
0.6 


Namur 


09 


3 


107 


33 


1 


42 


9 


1 33 4 


7.3 


1.8 




10 


4 


397 


30 


2 


295 


16 


1 80 12 


22.1 


1.7 


Walloon Brabant 


09 
10 
09 


1 


85 


51 










17.0 


10.2 


All (ten) provinces 


28 


4016 


576 


20 


3893 


268 


17 2019 259 


30.5 


3.4 




10 


31 


5761 


298 


23 


4068 


209 


6 968 65 


36.0 


1.9 



humidity and vapour pressure deficit on the sampling 
day. Vapour pressure deficit 3 and 7 days prior to sam- 
pling was also significant but at a lower level. The influ- 
ence of wind speed was only significant at the greatest 
temporal scale considered i.e. 7 previous days. The 



abundance of adults was significantly correlated only to 
wind speed averaged for 7 previous days. At site level 
(Table 3), both mean abundances of nymphs and adults 
were found significantly associated with the proportions 
of forest area and a set of landscape configuration 



2010 




% 
■84 
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0,4 
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I | Monthly mean sampled nymphs j^^j Monthly mean sampled adults 
Monthly mean temperature interval o Monthly mean relative humidity — □— Monthly mean vapour pressure deficit 



Figure 2 Temporal variation of ticks in 2009 and 2010. Red and blue lines indicate changes of monthly mean relative humidity and vapour 
pressure deficit. Grey bars indicate the interval between maximum and minimum temperatures. Mean values of ticks were calculated per 100 m 2 . 
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Table 2 Sampling level bivariate regressions 

Variables Descriptive Statistics Coefficient (Std. Error) in Regressions 







Mean 


Range 


Nymph abundance 


Adult abundance 


Daily minimum temperature (°C) 


Sampling day 


7.97 


-2.40-16.80 


0.03 (0.03) 


0.01 (0.02) 




3 previous days mean 


8.17 


0.03-18.07 


0.05 (0.03) 


0.03 (0.03) 




7 previous days mean 


8.21 


-1.20-18.56 


0.05 (0.03) 


0.03 (0.03) 


Daily maximum temperature (°C) 


Sampling day 


19.55 


8.40-33.00 


0.05 (0.03)* 


0.02 (0.02) 




3 previous days mean 


18.73 


7.40-29.80 


0.04 (0.03) 


0.01 (0.02) 




7 previous days mean 


18.44 


5.16-29.64 


0.03 (0.03) 


0.01 (0.02) 


Daily precipitation (mm) 


Sampling day 


1.28 


0-24.30 


-0.02 (0.04) 


-0.02 (0.03) 




3 previous days mean 


1.44 


0-21.47 


-0.02 (0.05) 


-0.04 (0.05) 




7 previous days mean 


1.93 


0-11.79 


0.04 (0.06) 


-0.01 (0.05) 


Daily wind speed (m/s) 


Sampling day 


3.13 


0.90-8.00 


-0.14 (0.10) 


-0.14 (0.09) 




3 previous days mean 


3.19 


0.90-7.60 


-0.12 (0.13) 


-0.10 (0.11) 




7 previous days mean 


3.27 


1.17-8.99 


-0.49 (0.1 4) *** 


-0.26 (0.12)* 


Daily relatively humidity (%) 


Sampling day 


71.75 


52.00-93.00 


-0.04 (0.01)* 


-0.01 (0.01) 




3 previous days mean 


74.36 


47.67-94.67 


-0.03 (0.02) 


-0.01 (0.01) 




7 previous days mean 


74.19 


51.00-91.86 


-0.03 (0.02) 


-0.001 (0.02) 


Daily vapour pressure deficit (kPa) 


Sampling day 


0.47 


0.11-1.25 


1.37 (0.55)* 


0.41 (0.48) 




3 previous days mean 


0.42 


0.07-1.35 


1.09 (0.62). 


0.18 (0.53) 




7 previous days mean 


0.42 


0.12-1.17 


1.19 (0.67). 


0.06 (0.58) 



Significant variables are highlighted in bold. 
***P < 0.001 ; **P < 0.01; *P < 0.05; .P < 0.1 . 



measures including forest edge density, shape index, 
nearest neighbour distance and aggregation index of for- 
est patches. As the buffer radius increased, the signifi- 
cance of all landscape composition and configuration 
characteristics increased. Site mean abundance of 
nymphs was also found negatively related to the percent- 
age of clay in the soil and the abundance categories of 
roe deer shot. 

Multi-level regressions 

Variables with P < 0.1 in single-level bivariate regressions 
were selected for multi-level regressions. A total number 
of 4 201 models for nymphs and 273 models for adults 
were generated by adding up to two variables for each 
level. Models containing significant collinearity and any 
non-significant covariates were removed, leaving 258 
models for nymph abundance and 60 models for adult 
abundance. Based on their improvement on model fit- 
ness (i.e. AAIC), there were 136 and 29 above-average 
models for nymphs and adults respectively. 

The relative importance of variables is presented in 
Table 4. The highest summed value of Akaike weights in 
above-average models of nymph abundance were 
reached by, at sampling level, daily relative humidity, 
mean value of daily wind speed in 7 previous days and 
daily maximum temperature and by, at site level, euclid- 
ean mean nearest-neighbour distance, area-weighted 



mean forest shape index and area-weighted mean patch 
fractal dimension. The most important predictors of 
adult abundance were the mean value of daily wind 
speed in 7 previous days at sampling level and area- 
weighted mean patch fractal dimension, forest edge 
density, euclidean mean nearest-neighbour distance and 
proportion of broad-leaved forest at site level. 

For all the top five nymph and adult models, neither 
spatial autocorrelation in residuals nor significant ran- 
dom slopes and cross-level effects were found. The 
between-site variances were greatly reduced in models 
with parameters than in the intercept-only model (model 
0, Tables 5 and 6). By including the site-level variables, 
site-level variances in the best models (model 1, Tables 5 
and 6) for nymph abundance and for adult abundance 
have respectively dropped by 75.9% and 84.5%. These 
findings indicate that the variability of nymph and adult 
tick abundance among sites can be adequately explained 
by the site-level variables included in the top models. 

Discussion 

Human risk of contracting tick-borne diseases is directly 
linked to the timing of the questing activity and the distri- 
bution of ticks. Therefore, it is essential to recognise rela- 
tions between abundance of questing ticks and factors 
associated with tick habitat and host pattern. In this study, 
we investigated these relations for Belgium using a multi- 
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Table 3 Site level bivariate regressions 



Vdl IdUltrb 




Descriptive Statistics 


Coefficient (Std. Error) in Regressions 






Mean 


Range 


Site mean nymph 
abundance 


Site mean adult 
abundance 


Landscape composition variables 












Proportion of broad -leaved forest (%) 


r a = 500 m 


24.98 


0-94.95 


2.84 (0.99)** 


1.67 (0.82)* 




r = 1 000 m 


21.66 


0.88-80.30 


4.14 (1.13)*** 


2.75 (0.93)** 




r= 1500 m 


21.27 


1 .06-85.32 


4.61 (1.14)*** 


3.62 (0.91)*** 


Proportion of coniferous forest (%) 


r = 500 m 


2.13 


0-29.18 


5.89 (5.01) 


3.67 (4.03) 




r = 1 000 m 


3.92 


0-86.68 


1.13 (1.97) 


0.63 (1.60) 




r= 1500 m 


2.50 


0-25.16 


14.61 (5.52)** 


10.25 (4.42)* 


Proportion of crops and pastures (%) 


r = 500 m 


35.25 


0-89.34 


-0.37 (1.00) 


-0.10 (0.82) 




r = 1 000 m 


35.72 


0.11-77.07 


-0.31 (1.25) 


-0.03 (1.02) 




r= 1500 m 


40.10 


6.70-83.70 


-0.61 (1.36) 


-0.90 (1.10) 


Proportion of heathland (%) 


r = 500 m 


0.11 


0-1.11 


1.78 (102.86) 


-16.94 (84.09) 




r = 1 000 m 


0.14 


0-1.19 


-31.32 (109.82) 


-23.92 (89.91) 




r= 1500 m 


0.34 


0-11.28 


-34.64 (16.62)* 


-18.81 (14.88) 


Landscape (forest) configuration variables 












The number of forest patches 


r = 500 m 


11.16 


0-27.00 


-0.03 (0.04) 


0.002 (0.03) 




r= 1000 m 


51.47 


6.00-104.00 


0.001 (0.01) 


0.005 (0.008) 




r= 1500 m 


124.90 


12.00-240.00 


0.001 (0.004) 


0.002 (0.003) 


Forest edge density (m/ha) 


r = 500 m 


21.56 


0-76.97 


0.03 (0.01)* 


0.02 (0.01). 




r = 1 000 m 


9.83 


0-48.32 


0.06 (0.02)*** 


0.06 (0.01)*** 




r= 1500 m 


14.12 


0-47.72 


0.07 (0.02)*** 


0.06 (0.01)*** 


Area-weighted mean forest shape index 


r = 500 m 


2.06 


0-3.5 


0.33 (0.37) 


0.42 (0.30) 




r= 1000 m 


2.67 


1.02-4.25 


0.66 (0.26)* 


0.79 (0.20)*** 




r= 1500 m 


3.01 


1.05-5.13 


0.72 (0.19)*** 


0.80 (0.14)*** 


Area-weighted mean forest patch fractal dimension 


r = 500 m 


1.11 


0-1.22 


0.01 (1.50) 


0.49 (1.24) 




r= 1000 m 


1.15 


1.02-1.23 


10.47 (4.32)* 


13.13 (3.45)*** 




r= 1500 m 


1.16 


1.03-1.24 


11.94 (3.77)** 


14.71 (2.99)*** 


Euclidean mean nearest-neighbour distance 
between forest patches (m) 


r = 500 m 
r = 1 000 m 


102.99 
89.03 


0-589.62 
63.88-339.93 


-0.02 (0.01)** 
-0.04 (0.01)** 


-0.005 (0.003). 
-0.08 (0.04)* 




r= 1500 m 


88.64 


63.81-195.46 


-0.04 (0.01)*** 


-0.03 (0.01)*** 


Aggregation index of forest patches 


r = 500 m 


63.13 


0-94.25 


0.02 (0.01) 


0.01 (0.01) 




r = 1 000 m 


59.42 


23.94-89.51 


0.04 (0.01)** 


0.04 (0.01)*** 




r= 1500 m 


57.67 


24.22-90.23 


0.04 (0.01)*** 


0.05 (0.01)*** 


So/7 variables 












Drainage degree 




2.98 


0-6 


-0.17 (0.17) 


0.06 (0.14) 


Proportion of gravel (%) 




3.92 


0-20.00 


-1.51 (3.08) 


-1.40 (2.52) 


Proportion of sand (%) 




37.35 


0-90.00 


1.01 (0.72) 


0.74 (0.59) 


Proportion of slit (%) 




41.92 


0-85.00 


-0.25 (0.86) 


-0.14 (0.70) 


Proportion of clay (%) 




10.92 


0-50.00 


-5.12 (2.45)* 


-3.18 (2.05) 


Wildlife variables 












Categories of shot roe deer number per 1 km 2 forest 




1.98 


1.00-5.00 


-0.31 (0.18). 


-0.20 (0.15) 


Number of shot roe deer per 1 km 2 forest 




12.09 


0-47.29 


-0.02 (0.01) 


-0.01(0.01) 



Significant variables are highlighted in bold. 

a r : Buffer zone radius. 

***P < 0.001 ; **P < 0.01 ; *P < 0.05; .P < 0.1 . 



Li et al. Parasites & Vectors 2012, 5:149 
http://www.parasitesandvectors.eom/content/5/1/149 



Page 8 of 1 1 



Table 4 Relative importance of variables (Sum of Akaike weights >0.1) in above-average models 


Level 


Variable 


Sum of Akaike weights 


Key determinants to nymph abundance: 






Sampling level 


Daily relative humidity (sampling day) 


0.56 




Daily wind speed (7 previous days mean) 


0.31 




Daily maximum temperature (sampling day) 


0.16 


Site level 


Euclidean mean nearest-neighbour distance (r a = 500 m) 


0.99 




Area-weighted mean forest shape index (r= 1500 m) 


0.70 




Area-weighted mean patch fractal dimension (r= 1500 m) 


0.24 


Key determinants to adult abundance: 






Sampling level 


Daily wind speed (7 previous days mean) 


0.77 


Site level 


Area-weighted mean patch fractal dimension (r= 1500 m) 


0.84 




Forest edge density (r= 1500 m) 


0.34 




Euclidean mean nearest-neighbour distance (r= 1500 m) 


0.22 




Proportion of broad-leaved forest (r= 1500 m) 


0.13 




Forest edge density (r= 1000 m) 


0.11 



r : Buffer zone radius. 



level approach. In general, /. ricinus ticks were found 
widely distributed in the forested areas of Belgium. The 
abundance of questing nymphs appeared sensitive to more 
selected variables than adults. The spatial pattern of ticks 
sampled is similar to that of recorded cases of Lyme dis- 
ease in Belgium, but not exactly the same [27]. Spatial var- 
iations of disease incidence depend not only on the 
distribution of the vector but also on the land use patterns 
that may influence human exposure to ticks [28]. 

Consistent with current knowledge, tick questing pat- 
terns in Belgium were influenced by short term weather 
conditions at sampling level. The dryness of the 



environment seemed related to different effects on the 
sampling results at daily and weekly scales. 

At the daily level, environment dryness was related to 
more sampled ticks. However, while relative humidity 
and vapour pressure deficit were highlighted for nymphs 
at the sample level, only relative humidity was retained, 
with a marginal predictive power, in the above-average 
models. This is in contrast to previous studies [2,3,9]. 
Similar to our findings, Hubalek et al. [29] and Ruiz- 
Fons and Gilbert [30] found a negative association be- 
tween daily abundance of questing nymphs and relative 
humidity. Tagliapietra et al. [14] also identified an 



Table 5 Null model and the top five models for nymph abundance 





Model 0 


Model 1 


Model 2 


Model 3 


Model 4 


Model 5 


Fixed effect 


Coefficient (Std. Error) 










Intercept 


3.29 (0.39)*** 


3.77 (1.14)*** 


3.94(1.38)** 


1.87 (1.05). 


-13.29 (5.35)* 


-16.28 (5.60)** 


Level- 1 : Sampling 














Daily wind speed (7 previous days mean) 




-0.33 (0.01)* 






-0.30 (0.15)* 




Daily relative humidity (sampling day) 






-0.02 (0.01). 




-0.02 (0.01). 


-0.02 (0.01). 


Daily maximum temperature (sampling day) 








0.03 (0.02). 






Level-2 : Site 














Euclidean mean nearest-neighbour distance (r a = 500 m) 




-0.02 (0.01)*** 


-0.02 (0.00)*** 


-0.02 (0.01)* 


-0.03 (0.01)*** 


-0.03 (0.01) ** 


Area-weighted mean forest shape index (r= 1500 m) 




1.01 (0.21)*** 


1.10 (0.23)*** 


1 .06 (0.22)* 






Area-weighted mean patch fractal dimension (r= 1500 m) 










18.63 (4.32)*** 


20.50 (4.65). 


Random effect 


Variance (Std. 


Error) 










Intercept 


5.46 (0.21) 


1.32 (0.10) 


2.00 (0.13) 


1.73 (0.12) 


1.62 (0.11) 


2.17 (0.13) 


AIC 


1433.25 


1391.81 


1393.49 


1393.92 


1394.68 


1396.54 


AAIC 


0 


41.44 


39.76 


39.33 


38.56 


36.71 



r : Buffer zone radius. 



; P < 0.001 ; **P < 0.01 ; *P < 0.05; .P < 0.1 . 
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Table 6 Null model and the top five models for adult abundance 





Model 0 Model 1 


Model 2 


Model 3 


Model 4 


Model 5 


Fixed effect 


Coefficient (Std. Error) 










Intercept 


1.65 (0.27)*** -12.95 (3.53)**^ 


* -12.95 (2.69)*** 


-13.61 (3.70)**^ 


* -9.62 (5.50). 


-17.08 (3.67)*** 


Level— 1 : Sampling 












Daily wind speed (7 previous days mean) 


-0.23 (0.14). 


-0.31 (0.12)** 


-0.27 (0.15). 


-0.23 (0.13). 




Level-2 : Site 












Proportion of broad-leaved forest (r a = 1 500 m) 




1.18 (0.44)** 






-0.24 (0.13). 


Forest edge density (r= 1000 m) 






0.02 (0.01). 






Forest edge density (r= 1500 m) 


0.03 (0.02)* 










Area-weighted mean patch fractal 
dimension (r= 1500 m) 


12.94 (3.04)* 


13.46 (2.32)*** 


13.71 (3.10)*** 


11.76 (4.07). 


16.85 (3.03) *** 


Euclidean mean nearest-neighbour 
distance (r= 1500 m) 








-0.02 (0.01)*** 




Random effect 


Variance (Std. Error) 










Intercept 


1.71 (0.12) 0.26 (0.05) 


1 .07e-07 (2.93e-5) 


0.27 (0.05) 


0.33 (0.05) 


0.47 (0.06) 


AIC 


881.85 853.93 


855.50 


855.78 


856.05 


856.76 


AAIC 


0 27.92 


26.35 


26.07 


25.80 


25.09 



r : Buffer zone radius. 



■P < 0.001 ; **P < 0.01 ; *P < 0.05; .P < 0.1 . 



optimal value of 0.98 kPa for vapour pressure deficit 
below which the effect was positive. While the reason 
remains unknown for such positive relation between 
short-term environment dryness and abundance of 
sampled nymphs, several explanations can be proposed. 
Firstly, lower vapour pressure deficit may lead to a 
longer period of questing [4] and a lower mortality rate 
[31]. This may contribute to a high feeding success of 
nymphs on their natural hosts and thus reduce the 
chance of encountering questing nymphs during sam- 
pling. Secondly, meteorological data derived from cli- 
mate stations may lead to uncertainty when used for 
forested sampling sites more than 10 km away. In a for- 
est ecosystem, the soil and ground-level humidity may 
be adjusted by ground vegetation or litter [32] and air 
temperature is modified by canopy [33]. These microcli- 
matic effects can respond to short-term variations of 
weather conditions. Hence, field measurement of me- 
teorological variables is required to gain more insight 
into the biology of ticks. Thirdly, humidity gradient 
within the vegetation layer may also influence the verti- 
cal distribution of ticks. The atmosphere is permanently 
moist near the roots where ticks remain quiescent 
(mostly shorter than 1 day [4]), but can be generally 
dryer at the vegetation tips where ticks quest. Ticks will 
need to rehydrate every now and then to avert desicca- 
tion and must return to the moist ground layer periodic- 
ally to do so. Ixodes ricinus nymphs are known to be 
able to establish "avoiding reactions" [34], i.e. when they 
go from the tip of vegetation toward the root and ap- 
proach the vertical humidity boundary, they may move 



to the dry side to avoid the change of humidity. This 
may keep ticks near the tip of vegetation till desiccation 
and thus increase the probability of being sampled. 

At the weekly level, drying of the environment seemed 
to have resulted in a lower abundance of sampled nymphs 
and adults. Wind speed in the previous 7 days was found 
to be important in a negative relationship to the daily 
questing abundance of both nymphs and adults. Similarly, 
in Texas, a negative relation was found between the aver- 
age wind speed of the 3 weeks prior to the sampling day 
and the abundance of immature Amblyomma cajennense 
ticks [35]. Continuous wind may increase the dispersal of 
carbon dioxide which is an attractant for Ixodid ticks [36] . 
Other than this, it may fasten the forest air replacement 
by dryer air from less vegetated areas. Such consequences 
on a relative long term (i.e. a week) can lead to an increase 
in tick mortality. 

In this study, sampling focused on favourable areas, 
forests, allowing comparison of specific features of 
forested landscape structure. Among identified site-level 
determinants, landscape composition factor (i.e. cover- 
age of broad-leaved forests) seemed to be important to 
the adult abundance only. However, forest fragmentation 
was important to the abundance of both nymphs and 
adults. In general, more ticks were sampled in land- 
scapes where forest is more fragmented, with longer 
edges and patches that are close to each other. In Bel- 
gium, as in many places, these features result from the 
way landscapes are shaped by human use of the land. In 
the field, edge density, shape index and fractal dimension 
are all good indicators to plant species richness [37]. 
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This may be attractive to foraging deer, and may there- 
fore contribute to high local tick abundance. Besides, 
smaller between-patch distance and lower aggregation 
index of forest can indicate a higher probability of hosts 
moving between the forest patches which, in theory, can 
contribute to a higher tick density in the forest [38,39]. 
Interestingly, nymphs and adults responded differently 
to the buffer radii used for the nearest neighbour dis- 
tance. It appeared to be more important for nymphs at a 
smaller geographical extent (buffer radius = 500 m, area 
= 79 ha) and for adults at a larger geographical extent 
(buffer radius = 1500 m, area = 707 ha). This may be 
associated to the home ranges of the preferred host for 
tick at different life stages. The distribution of nymphs 
can be influenced by the movements of hosts on which 
larvae feed, whereas the distribution of adults can be 
influenced by the movements of hosts on which nymphs 
feed. Larval ticks tend to feed on small-sized animals 
such as rodents whose home ranges are typically smaller 
than 79 ha. For example, the home ranges of bank vole 
and wood mice are around 0.2 ha [40,41]. Nymphs feed 
not only on rodents but also on larger animals such as 
roe deer and red deer whose home ranges can be as 
large as 100 ha and 750 ha respectively [42]. Such differ- 
ences can advise the design of the environmental data 
collection schemes for the future empirical works that 
specifically focus on tick life stages. 

Roe deer hunting was not identified as a key determin- 
ant to the dynamics of ticks from the multi-level models. 
The number of roe deer shot was found negatively but 
marginally associated with questing nymphal abundance 
in single level regressions only. Despite the fact that what 
proportion of the total population is shot is unknown, the 
hunting of roe deer seemed to reduce the abundance of 
nymphs. Hunting reduces both the roe deer population 
and the feeding success of adult ticks. This may result in a 
decrease of local tick population. Besides, hunting activity 
is performed mostly in relation to forests. Such forests 
normally contain various and abundant animal species 
that are suitable to host ticks. As a result, the natural tick 
feeding success can be high and only a few ticks would be 
left to be encountered during sampling. 

Several potential limitations of the present study 
should be noted. Firstly, though the abundance of 
sampled ticks is an important indicator of human expos- 
ure to questing ticks, it does not equal the tick density. 
The site-level sampled results thus cannot adequately re- 
flect the tick density differences in Belgium. Further in- 
vestigation on the relations between tick abundance to 
hosts populations and spatial dynamics including move- 
ments, would be useful. Secondly, the effects of wildlife 
management need to be further investigated. Hunting 
data of 2009 were used as a site-level variable. However, 
it is likely that there would be inter-annual variations in 



shooting patterns due to the change of wildlife and land 
management policy, or failed execution of local shooting 
plan. Also, the influence of such variation on tick popu- 
lation may appear only after a year or two [43]. Better 
data on wildlife, including large vertebrates, is therefore 
crucial for understanding better the links between tick 
abundance and host abundance. Thirdly, current inter- 
pretation of the sampled results can be modified by in- 
cluding microhabitat factors of the sampled area. Even if 
in the same forest stand, tick abundance has been shown 
to vary under different topographical, soil and ground 
vegetation conditions [6,15]. Understanding their effects 
can in turn contribute to a better standardisation of the 
sampling procedure in large-scale investigations. Finally, 
the small group sizes can significantly influence the pre- 
dictive power of multi-level models [22] and may result 
in an overestimation of site level variance [44]. Continu- 
ing studies for Belgium can be based on the present 
study and should preferably focus on establishing more 
seasonal samplings on fewer sites. 

Conclusions 

A multi-level approach was performed to study the spatio- 
temporal dynamics of L ricinus ticks in Belgium. Several 
environmental determinants were identified for daily abun- 
dance of questing ticks including temperature, relative 
humidity and wind speed at sampling level and forest frag- 
mentation at site level. The study showed that the effects of 
environmental variables on tick life stages responded differ- 
ently to different spatial and temporal scales. It can contrib- 
ute to optimise spatio-temporal data collection schemes for 
future empirical investigation and to introduce the potential 
parameters for the future biological tick models. 
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